I'm having a problem with a few different version of my run_mcmc function, where it's giving shitty contours. I'm gonna put it into a notebook so I can more easily make plots and step through different processes in it.


In [1]:
from pearce.emulator import OriginalRecipe, ExtraCrispy
from pearce.mocks.customHODModels import *
from pearce.mocks import cat_dict
import numpy as np
from os import path

In [2]:
from multiprocessing import cpu_count
import warnings
from itertools import izip

import numpy as np
import emcee as mc
from scipy.linalg import inv

In [3]:
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()

In [4]:
def lnprior(theta, param_names, *args):
    for p, t in izip(param_names, theta):
        low, high = _emu.get_param_bounds(p)
        if np.isnan(t) or t < low or t > high:
            return -np.inf
    return 0

In [5]:
def lnlike(theta, param_names, fixed_params, r_bin_centers, y, combined_inv_cov, obs_nd, obs_nd_err, nd_func_name):
    param_dict = dict(izip(param_names, theta))
    param_dict.update(fixed_params)
    y_bar = _emu.emulate_wrt_r(param_dict, r_bin_centers)[0]
    # should chi2 be calculated in log or linear?
    # answer: the user is responsible for taking the log before it comes here.
    delta = y_bar - y
    print 'y',y
    print 'ybar',y_bar
    #print y_bar
    chi2 = -0.5 * np.dot(delta, np.dot(combined_inv_cov, delta))

    return chi2# - 0.5 * ((obs_nd - getattr(_cat, nd_func_name)(param_dict)) / obs_nd_err) ** 2

In [6]:
def lnprob(theta, *args):
    lp = lnprior(theta, *args)
    if not np.isfinite(lp):
        return -np.inf

    return lp + lnlike(theta, *args)

In [7]:
def _run_tests(y, cov, r_bin_centers, param_names, fixed_params, ncores, nd_func_name):

    assert ncores == 'all' or ncores > 0
    if type(ncores) is not str:
        assert int(ncores) == ncores

    max_cores = cpu_count()
    if ncores == 'all':
        ncores = max_cores
    elif ncores > max_cores:
        warnings.warn('ncores invalid. Changing from %d to maximum %d.' % (ncores, max_cores))
        ncores = max_cores
        # else, we're good!

    assert y.shape[0] == cov.shape[0] and cov.shape[1] == cov.shape[0]
    assert y.shape[0] == r_bin_centers.shape[0]

    # check we've defined all necessary params
    assert _emu.emulator_ndim == len(fixed_params) + len(param_names) + 1  # for r
    tmp = param_names[:]
    assert not any([key in param_names for key in fixed_params])  # param names can't include the
    tmp.extend(fixed_params.keys())
    assert _emu.check_param_names(tmp, ignore=['r'])

    assert hasattr(_cat, nd_func_name)

    return ncores

In [8]:
def _resume_from_previous(resume_from_previous, nwalkers, num_params):
    # load a previous chain
    # TODO add error messages here
    old_chain = np.loadtxt(resume_from_previous)
    if len(old_chain.shape) == 2:
        c = old_chain.reshape((nwalkers, -1, num_params))
        pos0 = c[:, -1, :]
    else:  # 3
        pos0 = old_chain[:, -1, :]

    return pos0

In [9]:
def _random_initial_guess(param_names, nwalkers, num_params):
    """
    Create a random initial guess for the sampler. Creates a 3-sigma gaussian ball around the center of the prior space.
    :param param_names:
        The names of the parameters in the emulator
    :param nwalkers:
        Number of walkers to initiate. Must be the same as in resume_from_previous
    :param num_params:
        Number of params to initiate, must be the same as in resume_from_previous
    :return: pos0, the initial position of each walker for the chain.
    """

    pos0 = np.zeros((nwalkers, num_params))
    for idx, pname in enumerate(param_names):
        low, high = _emu.get_param_bounds(pname)
        pos0[:, idx] = np.random.randn(nwalkers) * (np.abs(high - low) / 6.0) + (low + high) / 2.0
        # TODO variable with of the initial guess

    return pos0

In [10]:
def run_mcmc(emu, cat, param_names, y, cov, r_bin_centers, obs_nd, obs_nd_err, nd_func_name, \
             fixed_params={}, resume_from_previous=None, nwalkers=1000, nsteps=100, nburn=20, ncores='all'):
    _emu = emu
    _cat = cat
    global _emu
    global _cat

    ncores= _run_tests(y, cov, r_bin_centers,param_names, fixed_params, ncores, nd_func_name)
    num_params = len(param_names)
    combined_inv_cov = inv(_emu.ycov + cov)

    sampler = mc.EnsembleSampler(nwalkers, num_params, lnprob,
                                 threads=ncores, args=(param_names, fixed_params, r_bin_centers, y, combined_inv_cov, \
                                                       obs_nd, obs_nd_err, nd_func_name))

    if resume_from_previous is not None:
        try:
            assert nburn == 0
        except AssertionError:
            raise AssertionError("Cannot resume from previous chain with nburn != 0. Please change! ")
        # load a previous chain
        pos0 = _resume_from_previous(resume_from_previous, nwalkers, num_params)
    else:
        pos0 = _random_initial_guess(param_names, nwalkers, num_params)
    return pos0, (param_names, fixed_params, r_bin_centers, y, combined_inv_cov, \
                                                       obs_nd, obs_nd_err, nd_func_name)

    sampler.run_mcmc(pos0, nsteps)

    chain = sampler.chain[:, nburn:, :].reshape((-1, num_params))


<ipython-input-10-9dbb0c2bea65>:4: SyntaxWarning: name '_emu' is assigned to before global declaration
  global _emu
<ipython-input-10-9dbb0c2bea65>:5: SyntaxWarning: name '_cat' is assigned to before global declaration
  global _cat

In [11]:
training_dir = '/u/ki/swmclau2/des/PearceLHC_wp_z_sham_emulator/'

em_method = 'gp'
split_method = 'random'

load_fixed_params = {'z':0.0}

In [12]:
emu = ExtraCrispy(training_dir,10, 2, split_method, method=em_method, fixed_params=load_fixed_params)

In [13]:
#Remember if training data is an LHC can't load a fixed set, do that after
fixed_params = {'f_c':1.0}#,'logM1': 13.8 }# 'z':0.0}

cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[1.0]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!

In [14]:
#mbc = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/mbc.npy')
#cen_hod = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/cen_hod.npy')
#sat_hod = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/sat_hod.npy')

#cat.load_model(1.0, HOD=(HSAssembiasTabulatedCens, HSAssembiasTabulatedSats),\
#                hod_kwargs = {'prim_haloprop_vals': mbc,
#                              'cen_hod_vals':cen_hod,
#                              'sat_hod_vals':sat_hod})
#cat.load_catalog(1.0)
cat.load(1.0, HOD='redMagic')

In [15]:
emulation_point = [('f_c', 0.2), ('logM0', 12.0), ('sigma_logM', 0.366),
                    ('alpha', 1.083),('logM1', 13.7), ('logMmin', 12.233)]
#emulation_point = [('mean_occupation_centrals_assembias_param1',0.6),\
#                    ('mean_occupation_satellites_assembias_param1',-0.7)]

em_params = dict(emulation_point)

em_params.update(fixed_params)
#del em_params['z']

In [16]:
#rp_bins =  np.logspace(-1.1,1.6,18) 
#rp_bins.pop(1)
#rp_bins = np.array(rp_bins)
rp_bins = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/rp_bins.npy')
rpoints = (rp_bins[1:]+rp_bins[:-1])/2.0

In [17]:
wp_vals = []
nds = []
for i in xrange(2):
    cat.populate(em_params)
    wp_vals.append(cat.calc_wp(rp_bins, 40))
    nds.append(cat.calc_number_density())

In [18]:
#y = np.mean(np.log10(np.array(wp_vals)),axis = 0 )
y = np.log10(np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/sham_vpeak_shuffled_wp.npy'))
# TODO need a way to get a measurement cov for the shams
cov = np.cov(np.log10(np.array(wp_vals).T))#/np.sqrt(50)

#obs_nd = np.mean(np.array(nds))
obs_nd = np.loadtxt('/nfs/slac/g/ki/ki18/des/swmclau2/AB_tests/sham_vpeak_shuffled_nd.npy')
obs_nd_err = np.std(np.array(nds))

In [19]:
param_names = [k for k in em_params.iterkeys() if k not in fixed_params]

nwalkers = 10
nsteps = 10
nburn = 0

In [20]:
pos0, args = run_mcmc(emu, cat, param_names, y, cov, rpoints,obs_nd, obs_nd_err,'calc_analytic_nd', fixed_params = fixed_params,\
        nwalkers = nwalkers, nsteps = nsteps, nburn = nburn)#,\
        #resume_from_previous = '/u/ki/swmclau2/des/PearceMCMC/100_walkers_1000_steps_chain_shuffled_sham_no_nd.npy')#, ncores = 1)

In [21]:
from corner import corner
corner(pos0, labels = param_names); plt.show()

In [22]:
for t in pos0:
    print lnlike(t, *args)
    print


y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 2.56963155  2.4547732   2.31934408  2.14747442  1.99309698  1.86687827
  1.7561729   1.6549791   1.55474039  1.43757182  1.30402759  1.17071475
  1.04612661  0.93573157  0.81512888  0.58880817  0.23433805]
-38.3018419197

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 2.71249502  2.58267256  2.43221085  2.24550714  2.07913669  1.94037773
  1.81927064  1.71187451  1.60318596  1.48623494  1.3722795   1.2644453
  1.15343637  1.03672015  0.89614207  0.67148658  0.29372361]
-59.5021256882

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 2.52177175  2.41034177  2.27798658  2.11087566  1.96905858  1.86548275
  1.77734792  1.68598673  1.58218657  1.46060099  1.32403554  1.18453062
  1.05393225  0.94206687  0.81957919  0.59600507  0.27420822]
-47.4652703033

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 3.58884308  3.48187494  3.36480535  3.18603648  2.96240609  2.70225478
  2.42399997  2.1671799   1.95908143  1.79936986  1.65683991  1.48783661
  1.28816639  1.09444582  0.90105306  0.60816509  0.14469855]
-30.9603916657

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 2.38699216  2.28324281  2.16827958  2.02419196  1.90051909  1.79808161
  1.70092698  1.60616203  1.50635377  1.39298752  1.27688542  1.16482946
  1.05278604  0.94654425  0.83740837  0.6366615   0.27508511]
-52.1835349166

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 2.31555192  2.20516457  2.10274524  1.97274009  1.87384603  1.8249844
  1.77670539  1.70300696  1.60540356  1.49602353  1.38299973  1.25509998
  1.10740093  0.9597248   0.78041502  0.53745316  0.22180955]
-30.5976381351

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 1.78980492  1.79302243  1.8461786   1.83303722  1.80004256  1.77716724
  1.73516804  1.65677239  1.55361053  1.44602434  1.3290614   1.19825035
  1.06540806  0.93896773  0.79056759  0.55134392  0.21683564]
-31.0433100187

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 3.07780981  2.8694145   2.64865412  2.40707084  2.20168094  2.0442939
  1.92361969  1.82512165  1.72389756  1.6130158   1.49573634  1.36648027
  1.23351155  1.10053993  0.94379555  0.71669129  0.43188544]
-105.778858766

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 2.36257036  2.24121259  2.13248749  1.99112934  1.87520256  1.80453219
  1.74944694  1.68479058  1.60147097  1.49949401  1.37664054  1.2507177
  1.12738911  0.99193774  0.82965081  0.61226407  0.369852  ]
-74.4367397345

y [ 2.80952448  2.65328823  2.48301516  2.30362902  2.14627957  1.95362842
  1.80430798  1.68158807  1.57263899  1.46153458  1.3476285   1.20096019
  1.05796409  0.88708418  0.66177206  0.39845162 -0.03907167]
ybar [ 3.6110835   3.48014945  3.35333711  3.16031305  2.91447941  2.6318688
  2.35048122  2.1152614   1.93053989  1.78597303  1.64063641  1.46031376
  1.25504145  1.05929945  0.87318662  0.59839015  0.16936722]
-30.2139839677

/u/ki/swmclau2/.conda/envs/hodemulator/lib/python2.7/site-packages/pearce/emulator/emu.py:445: UserWarning: One value for r is outside the bounds (0.097, 33.714) of the emulator.
  pname, plow, phigh))

In [23]:
print cov


[[  1.53112009e-04   1.40805711e-05   4.55467353e-05   8.45308972e-05
    5.72326603e-05   3.24262962e-05  -5.83476649e-05  -3.28622899e-05
   -6.59968407e-05  -8.59525275e-05  -3.09713914e-05  -4.64673545e-05
   -3.64807224e-05  -7.23893526e-05  -7.84130022e-05  -1.18995042e-04
   -3.70858347e-05]
 [  1.40805711e-05   1.29488525e-06   4.18859401e-06   7.77367703e-06
    5.26326150e-06   2.98200495e-06  -5.36580017e-06  -3.02210004e-06
   -6.06923790e-06  -7.90441378e-06  -2.84820819e-06  -4.27325651e-06
   -3.35486033e-06  -6.65710961e-06  -7.21105979e-06  -1.09430877e-05
   -3.41050800e-06]
 [  4.55467353e-05   4.18859401e-06   1.35489379e-05   2.51456853e-05
    1.70251886e-05   9.64595750e-06  -1.73568727e-05  -9.77565398e-06
   -1.96322983e-05  -2.55685824e-05  -9.21316216e-06  -1.38227975e-05
   -1.08520411e-05  -2.15339000e-05  -2.33257749e-05  -3.53978482e-05
   -1.10320458e-05]
 [  8.45308972e-05   7.77367703e-06   2.51456853e-05   4.66682701e-05
    3.15973134e-05   1.79020831e-05  -3.22128911e-05  -1.81427889e-05
   -3.64358890e-05  -4.74531313e-05  -1.70988515e-05  -2.56539457e-05
   -2.01404725e-05  -3.99651012e-05  -4.32906698e-05  -6.56954193e-05
   -2.04745461e-05]
 [  5.72326603e-05   5.26326150e-06   1.70251886e-05   3.15973134e-05
    2.13933409e-05   1.21208206e-05  -2.18101252e-05  -1.22837933e-05
   -2.46693567e-05  -3.21287131e-05  -1.15769830e-05  -1.73693124e-05
   -1.36363490e-05  -2.70588523e-05  -2.93104685e-05  -4.44798735e-05
   -1.38625376e-05]
 [  3.24262962e-05   2.98200495e-06   9.64595750e-06   1.79020831e-05
    1.21208206e-05   6.86729078e-06  -1.23569580e-05  -6.95962619e-06
   -1.39769122e-05  -1.82031582e-05  -6.55916881e-06  -9.84092766e-06
   -7.72594338e-06  -1.53307282e-05  -1.66064259e-05  -2.52009525e-05
   -7.85409500e-06]
 [ -5.83476649e-05  -5.36580017e-06  -1.73568727e-05  -3.22128911e-05
   -2.18101252e-05  -1.23569580e-05   2.22350293e-05   1.25231058e-05
    2.51499642e-05   3.27546436e-05   1.18025254e-05   1.77077007e-05
    1.39020119e-05   2.75860118e-05   2.98814940e-05   4.53464288e-05
    1.41326071e-05]
 [ -3.28622899e-05  -3.02210004e-06  -9.77565398e-06  -1.81427889e-05
   -1.22837933e-05  -6.95962619e-06   1.25231058e-05   7.05320311e-06
    1.41648413e-05   1.84479121e-05   6.64736130e-06   9.97324564e-06
    7.82982396e-06   1.55368603e-05   1.68297107e-05   2.55397965e-05
    7.95969866e-06]
 [ -6.59968407e-05  -6.06923790e-06  -1.96322983e-05  -3.64358890e-05
   -2.46693567e-05  -1.39769122e-05   2.51499642e-05   1.41648413e-05
    2.84470370e-05   3.70486633e-05   1.33497954e-05   2.00291187e-05
    1.57245173e-05   3.12024420e-05   3.37988539e-05   5.12911878e-05
    1.59853427e-05]
 [ -8.59525275e-05  -7.90441378e-06  -2.55685824e-05  -4.74531313e-05
   -3.21287131e-05  -1.82031582e-05   3.27546436e-05   1.84479121e-05
    3.70486633e-05   4.82511923e-05   1.73864179e-05   2.60853907e-05
    2.04791925e-05   4.06372294e-05   4.40187270e-05   6.68002768e-05
    2.08188845e-05]
 [ -3.09713914e-05  -2.84820819e-06  -9.21316216e-06  -1.70988515e-05
   -1.15769830e-05  -6.55916881e-06   1.18025254e-05   6.64736130e-06
    1.33497954e-05   1.73864179e-05   6.26487166e-06   9.39938438e-06
    7.37929534e-06   1.46428683e-05   1.58613279e-05   2.40702348e-05
    7.50169703e-06]
 [ -4.64673545e-05  -4.27325651e-06  -1.38227975e-05  -2.56539457e-05
   -1.73693124e-05  -9.84092766e-06   1.77077007e-05   9.97324564e-06
    2.00291187e-05   2.60853907e-05   9.39938438e-06   1.41021926e-05
    1.10713893e-05   2.19691567e-05   2.37972501e-05   3.61133318e-05
    1.12550325e-05]
 [ -3.64807224e-05  -3.35486033e-06  -1.08520411e-05  -2.01404725e-05
   -1.36363490e-05  -7.72594338e-06   1.39020119e-05   7.82982396e-06
    1.57245173e-05   2.04791925e-05   7.37929534e-06   1.10713893e-05
    8.69195773e-06   1.72476078e-05   1.86828126e-05   2.83519569e-05
    8.83613279e-06]
 [ -7.23893526e-05  -6.65710961e-06  -2.15339000e-05  -3.99651012e-05
   -2.70588523e-05  -1.53307282e-05   2.75860118e-05   1.55368603e-05
    3.12024420e-05   4.06372294e-05   1.46428683e-05   2.19691567e-05
    1.72476078e-05   3.42247379e-05   3.70726405e-05   5.62592972e-05
    1.75336970e-05]
 [ -7.84130022e-05  -7.21105979e-06  -2.33257749e-05  -4.32906698e-05
   -2.93104685e-05  -1.66064259e-05   2.98814940e-05   1.68297107e-05
    3.37988539e-05   4.40187270e-05   1.58613279e-05   2.37972501e-05
    1.86828126e-05   3.70726405e-05   4.01575223e-05   6.09407356e-05
    1.89927078e-05]
 [ -1.18995042e-04  -1.09430877e-05  -3.53978482e-05  -6.56954193e-05
   -4.44798735e-05  -2.52009525e-05   4.53464288e-05   2.55397965e-05
    5.12911878e-05   6.68002768e-05   2.40702348e-05   3.61133318e-05
    2.83519569e-05   5.62592972e-05   6.09407356e-05   9.24801392e-05
    2.88222359e-05]
 [ -3.70858347e-05  -3.41050800e-06  -1.10320458e-05  -2.04745461e-05
   -1.38625376e-05  -7.85409500e-06   1.41326071e-05   7.95969866e-06
    1.59853427e-05   2.08188845e-05   7.50169703e-06   1.12550325e-05
    8.83613279e-06   1.75336970e-05   1.89927078e-05   2.88222359e-05
    8.98269931e-06]]

In [24]:
np.diag(emu.ycov)


Out[24]:
array([  1.61832087e+02,   6.08813713e+01,   2.39505655e+01,
         9.39384523e+00,   3.83750197e+00,   1.58146552e+00,
         6.27034229e-01,   2.05122302e-01,   6.49350133e-02,
         2.93088017e-02,   1.64375795e-02,   8.45161588e-03,
         4.08487598e-03,   2.38226288e-03,   1.41069581e-03,
         9.64712579e-04,   1.13262294e-03])

In [ ]: